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Abstract 

2D materials are well-known to exhibit interesting phenomena due to quantum conhnement. Here, 
we show that quantum confinement, together with structural anisotropy, result in an electric-field- 
tunable Dirac cone in 2D black phosphorus. Using density functional theory calculations, we find 
that an electric field, Ugxt, applied normal to a 2D black phosphorus thin film, can reduce the direct 
band gap of few-layer black phosphorus, resulting in an insulator-to-metal transition at a critical 
field, Ec. Increasing Uext beyond Ec can induce a Dirac cone in the system, provided the black 
phosphorus film is sufficiently thin. The electric field strength can tune the position of the Dirac 
cone and the Dirac-Fermi velocities, the latter being similar in magnitude to that in graphene. We 
show that the Dirac cone arises from an anisotropic interaction term between the frontier orbitals 
that are spatially separated due to the applied field, on different halves of the 2D slab. When 
this interaction term becomes vanishingly small for thicker films, the Dirac cone can no longer be 
induced. Spin-orbit coupling can gap out the Dirac cone at certain electric fields; however, a further 
increase in field strength reduces the spin-orbit-induced gap, eventually resulting in a topological- 
insulator-to-Dirac-semi-metal transition. 


Introduction 

Two-dimensional (2D) layered materials, where inter¬ 
layer interactions are dominated by weak van der Waals 
(vdW) force^has attracted tremendous attention in na- 
noelectronicW Of particular interest are quantum con¬ 
finement effects in which the electronic structure of the 
2D layered material changes abruptly as a function of 
the number of layers. A well-known example is the in¬ 
direct to direct band g^ transition when 2D M 0 S 2 is 
reduced to a single layefl Here, we show that quantum 
confinement results in the emergence of an electrically 
tunable Dirac cone in a recently discovered 2D material 
- 2D black phosphorus. 

Black phosphorus, a layered semiconductor material, is 
the thermodynamically stable form of phosphorus. Dif¬ 
ferent experimental techniques, such as mechanicaP and 
liquicP exfoliation, have been employed to thin down the 
bulk to a monolayer. Experimentalists have succeeded 
also in measuring a thickness-dependent transport gam 
ranging from ^0.30 eV (bulk) to 1.0 eV (monolayer )P, 
and in fabricating field effect transistors based on ultra- 
thin black phosphoruP. Such transistors exhibit high 
mobility of ^ 1000 cm^/V • s and appreciably high on/off 
ratios, up to 10^, opening up a possible potential ap¬ 
plication in nanoelectronicP. In contrast to previously 
studied 2D materials such as graphene and M 0 S 2 , 2D 
black phosphorus layers are highly anisotropic ( see Fig¬ 
ure 1), with zigzag chains in one direction {y in the 
figure) and armchair chains in the other direction {x). 
This anisotropy adds another dimension of interest to 2D 
materials, resulting in direction-dependent optical and 
transport propertied 

In this work, we explore from first principles calcu¬ 
lations the effect of external electric fields applied per- 



FIG. 1: (Color online) The optimized structure of a bilayer 
(2L) black phosphorus, (a) and (b) show the top and side 
view of 2L black phosphorus. An external vertical electric 
field, Eext is applied along z-axis and 2U denotes the corre¬ 
sponding potential difference across the layers. A is the po¬ 
tential difference between the two vertically inequivalent sites 
in each layer. For the sake of simplicity we use | A| = U/2 for 
our tight-binding model of the bilayer system. In (a), hop¬ 
ping parameters of the tight-binding model are indicated, t 
and t' denote intralayer and interlayer hopping parameters 
respectively. In (b), dark and light yellow correspond to two 
different layers of black phosphorus. 


pendicular to thin films of black phosphorus. Similar to 
the Stark effect previously predicted for boron nitride 
nanotubes, two-dimensional transition metal dichalco- 
genides, and other ID and 2D material^^®^ we find that 
the electric field localizes the valence band maximum 
(VBM) and conduction band minimum (CBM) states 
at opposite surfaces of the slab. The resulting poten¬ 
tial difference between the VBM and CBM states results 
in a decrease in the band gap, eventually leading to an 
insulator-to-metal transition at a critical applied field, 
Ec. However, in contrast to other materials, we also find 
that for black phosphorus films below a certain critical 
thickness, a further increase in applied field strength re- 
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FIG. 2: (Color online) Electronic structure of 4L black phosphorus as a function of an external electric field, Eext- In (a) 
we present: (upper panel) the band structure and (lower panel) the charge density, respectively at the CBM (top) and VBM 
(bottom) for Egxt = 0. In (b) and (c), the same quantities are shown for Egxt = 0.31 and 0.72 V/A, respectively. Note that 
the critical held, Ec, is ~ 0.48 V/A for the 4L black phosphorus. In (c), where Eext > Ec, a Dirac cone appears along the ky 
direction, and a band gap of ~39 meV opens along the kx direction. Inset in (c): zoomed-in bandstructure around the F point 
near the Fermi level, Ep- Band inversion is evident in (c). 


suits in a highly anisotropic opening in the band gap that 
leads to a formation of a Dirac cone. This result is not ob¬ 
served for black phosphorus hlms that are too thick. We 
show that the emergent Dirac cone physics arises from 
anisotropic interaction terms between the electric-held- 
induced quantum-conhned VBM and CBM states in few 
layer black phosphorus. 


RESULTS AND DISCUSSION 

Electronic structure at zero Eext We begin our 
study by calculating the electronic structure of pris¬ 
tine black phosphorus at zero external electric field. 
Electronic structure calculations are performed by us¬ 
ing density functional theory (DFT) within the Perdew- 
Burke-Ernzerhof parametrization (PBE)f^ of the gener¬ 
alized gradient approximation (GCA) of the exchange 
and correlation potential, as implemented in the Quan¬ 
tum Espress^^ package. To treat van der Waals (vdW) 
interactions, we employ the vdW+DF approach for the 
exchange-correlation functionaP3(See Methods for fur¬ 
ther details). 

Bulk black phosphorus has an orthorhombic crystal 
structure; the optimized lattice parameters of bulk black 
phosphorus are a=3.35 A, 6=4.68 A, and c=11.42 A, 
in good agreement with the corresponding experimen¬ 


tal values (a=3.31 A, 6=4.37 A, and c=10.47 A)P^. Bulk 
black phosphorus is a semiconductor and we predict a 
direct band gap of 0.43 eV at the F point. Although 
this band gap value is slightly larger than the experi¬ 
mental value of 0.35 e\El, it is in good agreement with 
other vdW-corrected density functional theory calcula¬ 
tion For thin film black phosphorus, we find that the 
structural parameters do not change significantly. On the 
other hand, the band gap increases monotonically with 
decreasing thickness, reaching 0.99 eV for the monolayer. 
Interestingly, the band gap remains direct at the F point 
when the thickness reduces from bulk to monolayer limit. 

As previously mentioned, black phosphorus is a lay¬ 
ered compound consisting of puckered phosphorene layers 
with highly anisotropic bonding within each layer (Figure 
1). We denote multi-layered black phosphorus by nL-BP, 
where n is the number of layers. |^a) shows the electronic 
structure of 4L-BP, which exhibits a direct band gap of 
0.60 eV at the F point. The valence band maximum 
(VBM) and conduction band minimum (CBM) charge 
density are both delocalized over the layers. Our anal¬ 
ysis of the projected density of states indicates that the 
charge density at the VBM and the CBM originate from 
phosphorus 3s, 3px and 3pz orbitals (see Supplementary 
Figure SI). The Pz orbitals have the largest contribution 
(> 70%) to the VBM and CBM among all the orbitals. 

Electronic structure for Eext < Ec We now discuss 
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FIG. 3: (Color online) Variation of the (PBE+vdW) band gap 
with the applied electric field, Eext, for 2L (black circles), 3L 
(red squares), 4L (green triangles), 5L (blue stars), 6L (cyan 
triangles), 7L (magenta diamonds), and 8L (orange triangles) 
-BP. Inset: variation of the Stark effect coefficient, SnL, with 
number of layers. Note that we consider the linear region of 
the band gap versus Eext curve to compute SnL , and the fitted 
linear part is indicated by the corresponding dashed line. 


the response of the electronic structure of nL-BP to a 
static external electric field with magnitude F^ext, applied 
normal to the layers. As E’ext increases, the band gap de¬ 
creases monotonically while remaining direct at F, and 
eventually closes at a critical field, Ec (see Figure [^b) 
and Figure]^. This phenomenon essentially arises from 
the Stark effect: the applied electric field lifts the degen¬ 
eracy among distinct layers in the thin film. After the de¬ 
generacy is lifted, states that are localized in regions with 
higher potential will have a higher energy due to elec¬ 
trostatic interaction with the external field, while those 
localized in regions with lower potential will have a lower 
energy. Consequently, the VBM becomes localized on the 
surface layer with higher potential, while the CBM be¬ 
comes localized on the opposite surface layer with lower 
potential, resulting in a decrease in band gap. We note 
that, strictly speaking, the VBM and CBM states also 
have tails into the interior layers (Figure l^b-c)). These 
tails are important for the VBM and CBM states to inter¬ 
act when they become energy degenerate at the critical 
field (see later). 

Figure shows the modulation of band gap as a func¬ 
tion of F^ext for different nL-BP, where n runs from 2 to 
8. We find that the bandgap reduces very rapidly with 
E’ext as the thickness gets larger. Such a thickness depen¬ 
dence can be easily understood by considering the poten¬ 
tial drop across the layers. If one simply assumes that the 
potential drops linearly across the layers, a larger thick¬ 
ness requires a smaller electric field to sustain the same 
potential difference between the two edges. Furthermore, 
the band gap at zero field also decreases with increasing 
thickness. As such, the magnitude of the critical field, 
Ec, decreases with thickness. Indeed we observe that the 
value of Ec drops to 0.29 eV/A for 8L-BP. In practice. 


the potential drop is not strictly linear in space due to 
screening within the layers; however, an analysis of the 
screening effect shows that the potential drop in the in¬ 
terlayer region is much larger than that in the intralayer 
region, and the average potential in each layer drops lin¬ 
early across the layers (see Supplementary Figure S2). 

We quantify the above analysis as follows. Under 
the assumption that the potential drops linearly across 
the layers, the external field, F^ext, induces a poten¬ 
tial difference, AU, which can be expressed as AU = 
{eEc^^t/ k) {Z), where k is the dielectric constant along 
the direction of the applied field and (Z) is defined as 
(Z) = {Z)y — (^)q. Here {Z)y and {Z)^ denote the cen¬ 
ter positions (along the z-axis) of the localised VBM and 
CBM states, respective!}!!^. Since the VBM and CBM 
localize at opposite surfaces for large F^ext? {Z) = nc—2S, 
where c is the layer-layer distance and 6 is the position of 
(Z)y Q from the surfaced. We note that we are assum¬ 
ing here that S is independent of the applied field and 
film thickness. Although this is not likely to be valid, 
deviations should be small compared to (Z). Under this 
model, we can write the change in band gap, AF^g, as 
follows: 


AEg = eAV = (Z). (1) 

In general for the critical field, E^ we obtain from equa¬ 
tion [T] that 

-E^ = —{Z) = —{nc-2S). (2) 

K Kj 

Indeed, we find that EgjEc versus n is approximately 
linear (see Supplementary Figure S3), with deviations 
resulting from the assumption that S is constant. 

In the initial response to F^ext? the band gap varies 
non-linearly (likely quadratic) with U’ext, and the slope 
of the curve vanishes at F^ext = 0 (Figure 3). This is 
consistent with the fact that the variation in band gap 
should be independent of the polarity of the applied field, 
so that in this low field regime, the above linear analysis 
breaks down. However, for larger electric fields, for in¬ 
stance (F^ext/^c)> 0-1 for a specific thickness, AEg is an 
approximately linear function of F^ext (Figure]^ and this 
linear slope can be described as 


dEcy^l 


— eF'nL, 


( 3 ) 


where e is the electron charge and is defined as the 
linear giant stark effect (GSE) coefficient. Combining 
Equation [^with Equation the GSE coefficient can be 
written as 


= <^ = ( 4 ) 

K, 

Indeed, the inset of Figure shows that the Stark coef¬ 
ficient increases approximately linearly as a function of 
the number of layers. 
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FIG. 4: (Color online) (a) The PBE+vdW bandstructure of 4L-BP is plotted in the 2D brillouin zone for Ee^t = 0.72 V/A, 
above the critical field. The Dirac cone appears at ±A. (b) The PBE+vdW bandstructure along Y-P-X, for the three different 
Eext above the critical value, namely Eext = 0.51 (black), 0.72 (red), and 0.93 (green) V/A. (c) The bandgap that opens along 
T-X direction, Eg~^ , is plotted against (Eext — Ec) for 4L-BP (blue line) and 8L-BP (red line). Note that the values of Ec 
for 4L-BP and 8L-BP are 0.48 and 0.28 V/A, respectively, (d) The Eermi velocity along kx, Vp (red dashed line), and along 
ky, Vp (red solid line) are plotted as a function of (Eext — Ec) for 4L-BP. r'p is calculated as = SE/Sk at the Dirac point. 


At the critical field, the system becomes metallic when 
the VBM and CBM touch at E. Unlike the electrically 
driven Stoner magnetism previously reported for M 0 S 2 
nanoribbonJ^, there is no spin polarization in this system 
(Supplementary Figure S4). This is because the density 
of states at the Fermi level is much smaller in the 2D 
system than in the ID nanoribbons. 

Electronic structure beyond critical field An in¬ 
teresting and unique effect is observed in very thin-layer 
BPs (for example in 4L-BP), when an electric field larger 
than the critical field is applied. Specifically, a Dirac cone 
appears at a /c-point along F-Y, which we shall call A, as 
shown in Figure [^c) and Figure |^a-b). The valence 
band manifold and conduction band manifold touch at 
exactly two points, ±A, in the 2D Brillouin Zone [Fig¬ 
ure |^a)]. By plotting the charge density of the highest 
occupied and lowest unoccupied states at k-points close 
to the Dirac point (see Supplementary Figure S5), we see 
that the bands are inverted at the Dirac point, indicating 
that we have a new topological phase. As black phospho¬ 
rus films exhibit a direct band gap even in the presence 
of the external electric fields, the band inversion in itself 
may not be too surprising. However, the fact that the va¬ 
lence and conduction band manifolds touch at only one 
point, the Dirac point, is intriguing and deserves further 
analysis. 

As F^ext increases beyond the critical value, the Dirac 
point, A, shifts progressively from F toward Y in the 2D 
Brillouin zone [Figure Qb)]. This shifts occurs simulta¬ 
neously with an enhancement of the corresponding Dirac 
Fermi velocity, v-p [Figure [^c)]. The existence of the 
Dirac point is associated with an asymmetry between 
the F-Y and F-X directions - so that the valence and 
conduction bands touch only at the Dirac points and 
not in a ring in the 2D Brillouin Zone; neither does a 
band gap open up throughout the 2D Brillouin Zone. 
We quantify this asymmetry and therefore also the ro¬ 
bustness of the Dirac electronic structure by the band 


gap along F-X, E^~^. As shown in Figure Efc), 
increases initially as E’ext increases beyond me critical 
field; however, as F^ext continues to increase, de¬ 

creases. For 4L-BP, E^~^ reaches a maximum value of 
^39 meV at £’ext=0.72 V/A , where the corresponding 
X and ^-components of vp are 0.84x 10^ and 0.56x10^ 
m/s, respectively. These Fermi velocities are as large as 
the predicted value of ^ 0.86 x 10^ m/s in graphen^^. 
Interestingly, we find that the maximum value of E^~^ 
reduces to 8 meV when the thickness of black phospho¬ 
rus increases to 8L [FigureQc)], and vanishes for the 12L 
case (see Supplementary Figure S6). Salient features that 
arise from the above analysis are that F^g cannot be 
arbitrarily increased with electric field strength, and that 
E’g is smaller for thicker black phosphorus films. To 
understand the above intriguing phenomena, we develop 
a simple tight-binding model that captures the essential 
physics of black phosphorus thin films in the presence of 
an external electric field. 

Tight-binding model For simplicity we consider a 
tight-binding model for BP films that are two-layer thick 
(bilayer BP). Since we are concerned with the valence 
band (VB) and conduction band (CB) which stem from 
Pz orbitals, we include only Pz orbitals in our model. In 
the presence of large external electric fields, the VB and 
CB are mainly localized at the opposite surface layers 
of the thin film; therefore, the bilayer model can effec¬ 
tively also represent the VB and CB electronic structure 
in thicker films. 

We begin with a standard tight-binding Hamiltonian: 

H = Y^{ei + Ui) Cic\ + Uj Cic] + Y *ij (5) 

i <ij> <bi> 

where Ci is the energy of i-th site, Ui is the potential 
shift at the i-th site upon the applied vertical field (see 
figure [^, cl(cj) is the creation (annihilation) operator 
of electrons at site i{j), and tij and tC are respectively 
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FIG. 5: (Color online) The bandstructure of bilayer black 
phosphorus calculated by using a simple tight-binding model. 
In (a) the bandstructure is plotted for F^ext = 0 (black line) 
and F^ext <Ec (red line). In (b) the bandstructures are plotted 
for three different F^ext above the critical value, where Ui<U 2 
< 1/3 < 1/4 ■ For the applied bias U2 (green lines) and U3 (ma¬ 
genta lines), the interlayer hopping, t', is finite and kept fixed, 
whereas t' = 0 for the applied bias U 4 (cyan dashed lines). 
The horizontal dotted line corresponds to the Fermi level, 
(c) The kx — ky plot of the magnitude of < 0v|iF(tOl^c >, 
where 0v and 0c are the eigenvectors of the Hamiltonian H 
at = 0 for the highest occupied and the lowest unoccupied 
energy eigenvalues, respectively. We use eV for the 2D 

plot. 


the intralayer and interlayer inter-site hopping integrals 
as depicted in figure In the first instance, we include 
hopping terms up to the fourth and fifth nearest neighbor 
for the interlayer and intralayer interactions respectively, 
using the parameters reported previously in Ref.l^. This 
provides an accurate description of the VB and CB, as 
given by the GW approximation. Figure [^a) shows that 
our model nicely reproduces the band dispersion of the 
VB and CB in comparison with that of the GW bandJ^. 

When an external field is applied, the on-site energy 
of each site, 5^, is rigidly shifted by the corresponding 
potential difference, Ui (see Figure [^. For small Ui [ 
Figure [^a)] we see that the band gap closes as described 
above. When Ui increases above a critical value, Vc, 
the Dirac cone emerges along the F-Y direction, just as 



FIG. 6: (Golor online) Effect of spin-orbit coupling on the 
band structure of 4L black phosphorus, zoomed-in at the A 
point. Black and red lines denote the bandstructure calcu¬ 
lated without and with spin-orbit coupling, respectively, at 
F/ext — (a) 0.56, (b) 0.61 and (c) 0.66 V/A. 


was predicted by DFT, with the position of the Dirac 
cone shifting with the size of Uc (Figure |^b), green and 
magenta curves). 

However, interestingly, when the interlayer interaction 
terms t' are set to zero, U > Uc causes the band gap to 
vanish all over the two-dimensional Brillouin Zone (Fig¬ 
ure ib), cyan curve). This finding indicates that the in¬ 
terlayer interaction is critical to the emergence of a Dirac 
cone. Intuitively, we expect the VB and CB states to 
interact via these interlayer terms, since the states are 
localized at opposite surfaces of the thin film. The open¬ 
ing of a band gap everywhere in the Brillouin Zone ex¬ 
cept along F-Y thus results from interactions between the 
VB and CB states, mediated by the interlayer hopping 
integrals. This also explains why the maximum value 
of decreases as the thickness of the BP film in¬ 

creases and the interactions between the VB and CB be¬ 
come smaller. Thus, while the critical field Ec may be¬ 
come smaller with thicker films due to a decrease in the 
band gap with increasing thickness, the Dirac cone phe¬ 
nomenon will be absent beyond a critical thickness even 
with gap closure. We find that within DFT, when the 
film thickness increases to 12L, becomes negligible 

(equal to the precision of our calculations) (see Supple¬ 
mentary Figure S6). Quantum confinement is therefore 
important for the Dirac cone phenomenon. 

The unique Dirac cone band structure implies also that 
for U > Uc^ there is no interaction between the VB and 
CB states along the F-Y direction, but non-zero interac¬ 
tions along all other directions. This is a direct conse¬ 
quence of the anisotropic interlayer interactions in black 
phosphorus. Indeed, as shown in Figure l^c), the k- 
dependent Hamiltonian overlap matrix elements between 
the VB and CB states are highly anisotropic, going to 
zero along = 0, i.e. along the F-Y direction, but non¬ 
zero elsewhere in the Brillouin Zone. 

From Figure l^c), we can also see that for a fixed elec¬ 
tric field strength, the interaction between the VB and 
CB states increases as we move along F-X away from 
the F point. As the electric field strength is increased, 
the VB and CB cross further away from F. This ex- 
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FIG. 7: (Color online) Total energy per unit cell plotted 
against d — do for different values of F^ext in V/A. d is the 
inter-layer distance and do=3.59 Ais the equilibrium inter¬ 
layer distance. Inset shows the variation of binding-energy 
(BE), calculated within the PBE+vdW functional, as a func¬ 
tion of the applied electric field, Eext, for 2L black phosphorus. 



CO (eV) 


plains why initially increases as F^ext is increased [ 

Figurej^c)]. However, a larger electric field strength also 
implies greater localization of the VB and CB states, and 
therefore a smaller overall interaction term. These two 
effects compete with each other, resulting in an initial in¬ 
crease and then a subsequent decrease of as F^ext 

increases [ Figure |^c)]. 

Thus, we have shown that the Dirac cone phenomenon 
arising for F^ext > stems from a combination of quan¬ 
tum confinement and anisotropy effects in black phos¬ 
phorus. 

Bandstructure with spin-orbit correction As al¬ 
ready mentioned before, the band gap remains direct at 
the F point until the the gap closes completely. As a con¬ 
sequence, the band gap does not change significantly by 
the inclusion of spin-orbit coupling (SOC). However, the 
role of SOC becomes preeminent when E’ext > ^c- Inter¬ 
estingly, Figure ^a) shows that the Dirac cone is gapped 
out in the band structure by including the SOC, which 
demonstrates that the system is likely to be a quantum 
spin Hall insulator. The prediction can be explicitly con¬ 
firmed by calculating the Z 2 topological invariant of the 
systenPSl^ as shown in Ref For 4L-BP, the maxi¬ 
mum SOC-induced gap-opening at the Dirac point is as 
large as 3 meV at F^ext= 0.56 V/A. When F^ext is further 
increased, the SOC-induced gap decreases (to 2meV at 
^ext= 0.61 V/A) and subsequently, the Dirac semimetal 
phase appears at F^ext=0.66 V/A along the F-Y direction 
(Figurej^. Thus, the external electric field can be used to 
tune the opening and closure of the gap along F-Y when 
SOC is included. It is interesting to mention that the 
band inversion is controlled by electrically driven band 
overlap between the VB and CB, unlike other topologi¬ 
cal insulators, for example Bi 2 Se 3 , where band inversion 
is caused by spin-orbit coupling. 


FIG. 8: (Golor online) Tuning of optical properties by thick¬ 
ness and applied Eext- (a) The optical conductivity of 2L 
(black), 4L (red), 6L (green), and 8L (blue) BP, plotted as 
a function of frequency. Solid and dashed lines denote two 
different components of optical conductivity, axx and ayy , re¬ 
spectively. (b) axx (solid lines) and ayy (dashed lines) of 4L- 
BP are plotted against the frequency for two different applied 
fields, Eext =0.0 (black) and 0.21 (red) V/A. 


Robustness of results against functional Our re¬ 
sults above are robust against the exact choice of func¬ 
tional. Specifically, with the PBE functional, the band 
gap also decreases monotonically with thickness (see sup¬ 
plementary figure S7), but the rate of reduction of the 
PBE-bandgap is different from that of the vdW-corrected 
one. However, the qualitative physics, i.e. both the 
thickness-dependent bandgap, and the external field in¬ 
duced bandgap-closure and Dirac cone, consistently ap¬ 
pear regardless of the functional. Notably, the strength 
of critical field required for the insulator-to-metal tran¬ 
sition of bilayer BP is almost the same for both PBE 
and PBE-fvdW functionals. Moreover, although many- 
electron GW calculations are in principle required for 
accurate calculation of the band gap, the PBE and vdW 
corrected functionals underestimate and overestimate the 
bulk band gap respectively, and in both cases, the same 
physics emerges. 

Binding energy We also calculate the total energy as 
a function of interlayer distance in bilayer BP for differ¬ 
ent F^ext, as shown in Figure We find that the bind¬ 
ing energy reduces by increasing the strength of F^ext and 
eventually vanishes for a field of F^ext ^ 1.6 V/A, at which 
point the layers becomes unbounded. This field is much 
larger than the critical fields at which the Dirac cone 
phenomenon was predicted. Notably the equilibrium dis¬ 
tance between the layers is almost independent of F^ext- 
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This result suggests that an external field can trigger 
the exfoliation process of black phosphorus as an alter¬ 
native to mechanicaP and chemical exfoliatiorP. Note 
that few-layer graphene has been successfully exfoliated 
in the presence of ^ext^- 

Optical Properties Given the predicted changes in 
band gap as a function of thickness and electric field 
strength, we also consider a related quantity - the optical 
conductivity of the system. Figure [^a) shows both the 
XX and ^^-components of the optical conductivity, 
and <7yy^ as a function of frequency for 2L, 4L, 6L and 8L 
structures. We find that the onset of cjxx and cFyy depend 
on the dipole transitions at a point near R, and at the F 
point, respectively. As the bandgap reduces with increas¬ 
ing thickness, the onset of (Tyy also decreases, while that 
of (Jxx decreases slightly. By applying an external field to 
4L-BP, we find that the electric-field induced reduction 
in direct bandgap decreases the onset of cr^^^hereas that 
of cFxx remains almost unchanged [ FigureWb)]. 


CONCLUSIONS 

In summary, we have shown using first principles cal¬ 
culations that the application of an external electric field 
to few-layer BP results in a reduction in the direct band 
gap of the material, eventually resulting in an insulator- 
to-metal transition at a critical field. Above this critical 
field, a Dirac cone emerges at A along the T — Y direc¬ 
tion, provided that the BP film is sufficiently thin. We 
quantify the robustness of the Dirac cone by the energy 
gap in the T — X direction, • Using a tight-binding 

model, we show that quantum confinement and struc¬ 
tural anisotropy are critical for the emergent Dirac cone 
phenomenon. These conclusions further explain the sub¬ 
tle dependence of on the electric field strength and 

film thickness, providing us with a handle to tune the 
resulting electronic structure. The electric field strength 
can be used to tune the position of A and the Dirac-Fermi 
velocities, the latter being similar in magnitude to those 
in graphene. We also show that by including spin-orbit 
coupling, the Dirac point can be gapped out at certain 
electric fields, but increasing the electric field strength 
further can change the material from a topological insu¬ 
lator to a Dirac semi-metal. We have thus shown that the 
electronic band structure of few-layer black phosphorus 
can be tuned in very interesting ways by an applied elec¬ 
tric field. These results will motivate the search for other 
direct band gap materials where quantum confinement 
and structural anisotropy may lead to tunable emergent 
Dirac cone physics. 


METHODS 

Electronic structure calculations are performed by us¬ 
ing density functional theory within the generalized gra¬ 
dient approximation (GGA) of the exchange and cor¬ 


relation potential Perdew-Burke-Ernzerhof parametriza- 
tion (PBE]P^as implemented in the Quantum Espress^^ 
package. To treat van der Waals (vdW) interactions, 
we employ the vdW+DE approach for the exchange- 
correlation functionaP^. The wavefunction and the 
charge density are expanded using energy cutoffs of 50 
and 300 Ry, respectively. Brillouin zone sampling is done 
by using a (10x14x1) Monkhorst-Pack k-grid. Periodic 
boundary conditions have been included and a vacuum 
layer of at least 20 A is included to suppress the inter¬ 
action between the periodic images. The conjugate gra¬ 
dient algorithm is used to obtain optimized geometries, 
where all the atoms in the unit cell are allowed to relax 
until the forces on each atom are less than 0.002 eV/A, 
and the coordinates are fixed in the calculations using 
the finite external electric field. We have checked that 
for 4L-BP, the results are qualitatively similar when the 
atoms are relaxed in the presence of the applied electric 
field (see Supplementary Eigure S8). 

Simulations with an external electric field are per¬ 
formed using a periodic sawtooth-type potential perpen¬ 
dicular to the layer^^. In order to remove the effect of 
electrostatic interactions between the periodic supercells, 
a dipole correction scheme suggested by L. BengtssoiP^, 
is included. 

The real part of the optical conductivity is computed, 
using the following formula: 5R(cr) = uj x ^(e). Here the 
imaginary part of the dielectric tensor, can 

be expressed as 


^(Ca/3(^)) = 


47re^ 




-A[q//3 X f 
BZ ^ k{Ev — Ec) 


6 (^E\^ q E\^ y T Huj^ T 6 (^E\^^q E\^ y 


(6) 


where e and m are the electronic charge and mass, re¬ 
spectively, D is the lattice volume, W is the total num¬ 
ber of k-points, the indices c and v denote the conduction 
and the valance band, respectively, and / is the Eermi- 
Dirac distribution. is defined as follows: = 

(Mk,c|Pa|'Wk,«)(Mk,«|p^|Mk,c), where |Mk) and P are the 
single particle Bloch function and the dipole moment op- 
erator, respectively. An uniform k mesh of 36x26x1 is 
used for self consistent steps in order to calculate the 
optical conductivity. Note that our calculation does not 
include many-body effects which are required to obtain 
an accurate description of the optical conductivity. 
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Projected density of states (PDOS) 
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FIG. 9: Supporting Figure SI: The band-decomposed projected density of states (PDOS) on different orbitals for 4 layer 
black phosphorus. The densities of states are projected on 3s, 3px, Spy and 3pz orbitals. Upper, middle and lower panels show 
the PDOS at Ugxt = 0.00, 0.21, and 0.72 V/A, respectively. Note that the critical field is Ec = 0.48 V/A. At Uext = 0, the 
PDOS is summed over all the atoms in the unit-cell. For finite fields, for visualization purposes, we use the negative colour 
scale for PDOS that comes from the two layers kept at positive potential, and the positive colour scale for PDOS coming from 
the other two layers kept at negative potential, as indicated in Figure 1 in the main text. 
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FIG. 10: Supporting Figure S2: (a) Electric field induced charge density distribution across the layers, Ap = p(Eext) — 
p(Eext = 0), for 4L-BP and different values of Eext- p(Eext) is the charge density for an external held, Eext- p(Eext) is averaged 
over the transverse direction (xy-plane) and plotted along the out-of-plane direction (z-axis). The dashed vertical lines denote 
the positions of each atomic plane in 4L-BP. (b) Effective electrostatic potential across the layers, Af/gff = U{Ee^t) — U{Ee^t = 0), 
for 4L-BP and different values of Eext- 
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FIG. 11: Supporting Figure S3: (a) The polarization (dipole moment per unit area) , P, is plotted against the applied 
external field strength for 2L (black), 3L (red), 4L (green), 5L (blue), 6L (yellow), 7L (brown) and 8L (grey)-BP. (b) {EgjEc 
is plotted as a function of the number of layers. 
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(a) 



(b) DOS (states/eV) 



FIG. 12: Supporting Figure S4: (a) Spin-polarized bandstructure and (b) density of states of 4L-BP. The blue dashed line 
denotes the Fermi level. Red and black lines denote the spin-up and spin-down component of the DOS, respectively. 
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FIG. 13: Supporting Figure S5: The local density of states (LDOS) of highest occupied states (HOS) left) and lowest 
unoccupied states (LUS) (right) at the different k-points around the Dirac point, A, as indicated in the bandstucture (top). 
The figure shows that band inversion occurs at the A point. 
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FIG. 15: Supporting Figure S7: Comparison of electronic structure for different functionals, (a) The variation of the band- 
gaps calculated within PBE(black line) and vdW-corrected(red line) exchange correlation functionals, plotted as a function of 
the number of layers. Note that the PBE band gap of bulk BP is 0.12 eV, while the vdW-corrected band gap is 0.34 eV. The 
dashed line denotes the value of the experimental band gap (0.34 eV) of bulk BP. (b) The variation of the band gap of bilayer 
BP, calculated within PBE(black line) and vdW-corrected(red line) functionals, plotted against the applied field strength. 
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FIG. 16: Supporting Figure S8: The bandgap of 4L-BP as a functional of F^ext using the unrelaxed (black) and relaxed 
(red) coordinates under an applied field F^ext- 







